/*
 * Main.cpp
 *
 *  Created on: 13 Jul 2011
 *      Author: allan
 */

#include <iostream>
#include <iomanip>
using namespace std;

#include <Eigen/Core>

#include <GeoReact.h>

#include "KineticReactionPath/Example_CalciteDissolution.h"
#include "KineticReactionPath/Example_CalciteMagnesiteDolomite.h"
#include "KineticReactionPath/Example_CalciteCO2Dissolution.h"
#include "Solubility/Example_CO2Solubility.h"
#include "TRISpH/TRISpH.h"
#include "Test/Test_Compare_EquilibriumSolver_TotalEquilibriumSolver.h"

int main()
{
	//TRISpH(25.0, 1.0);
	
	//Example_CO2Solubility();
	//Example_CalciteCO2Dissolution();
	//Example_CalciteDissolution();
	//Example_CalciteMagnesiteDolomiteEquilibrium();
	
	Test_Compare_EquilibriumSolver_TotalEquilibriumSolver();
}

//class Problem : public NewtonProblem
//{
//public:
//	void Function(const VectorXd& x, VectorXd& F)
//	{
//		F[0] = x[0]*x[0] + x[1]*x[1] - 13;
//		F[1] = x[0]*x[1] + 6;
//	}
//	
//	void Jacobian(const VectorXd& x, MatrixXd& J)
//	{
//		J(0, 0) = 2*x[0];
//		J(0, 1) = 2*x[1];
//		J(1, 0) = x[1];
//		J(1, 1) = x[0];
//	}
//};
//
//int main()
//{
//	Problem problem;
//	
//	VectorXd x(2), F(2); x << 1, 4;
//	
//	NewtonSolver solver(2);
//	
//	solver.Solve(problem, x);
//	
//	problem.Function(x, F);
//	
//	cout << "x = " << x << endl;
//	cout << "F = " << F << endl;
//	
//	Example_CalciteMagnesiteDolomite();
//	//Example_CalciteMagnesiteDolomiteEquilibrium();
//}

//int main()
//{
//	GeoManager manager;
//	
//	manager.SpeciateAqueousSpecies({"H2O", "Calcite"});
//	manager.SetGaseousSpecies({});
//	manager.SetMineralSpecies({"Calcite"});
//	
//	GeoSystem system = manager.AssembleGeoSystem();
//	
//	cout << system << endl;
//	
//	GeoEquilibrator equilibrator(system);
//	
//	equilibrator.ImposeChargeBalance();
//	equilibrator.ImposeSpeciesMass("H2O", 55.508);
//	equilibrator.ImposeSpeciesMass("Ca[2+]", 1.0E-7);
//	equilibrator.ImposeSpeciesMass("HCO3[-]", 1.0E-7);
//	
////	equilibrator.ImposeChargeBalance();
////	equilibrator.ImposeAcidity(6.0);
////	equilibrator.ImposeSpeciesMass("H2O", 55.508);
////	equilibrator.ImposeElementMass("Ca", 2.0);
//	
//	GeoState state(system);
//	
//	equilibrator.Equilibrate(state);
//	
//	cout << state << endl;
//}

//int main()
//{
//	cout << ElementAtoms("H2O", "H") << endl;
//	cout << ElementAtoms("H2O", "O") << endl;
//	cout << ElementAtoms("Ru(OH)2SO4(aq)", "Ru") << endl;
//	cout << ElementAtoms("Ru(OH)2SO4(aq)", "O") << endl;
//	cout << ElementAtoms("Ru(OH)2SO4(aq)", "H") << endl;
//	cout << ElementAtoms("Ru(OH)2SO4(aq)", "S") << endl;
//	cout << ElementAtoms("Fe3(OH)2SO4(aq)", "Fe") << endl;
//	cout << ElementAtoms("Fe3(OH)2SO4(aq)", "F") << endl;
//	cout << endl;
//	cout << ElementAtoms("(UO2)11(CO3)6(OH)12[2-]", "U") << endl;
//	cout << ElementAtoms("(UO2)11(CO3)6(OH)12[2-]", "O") << endl;
//	cout << ElementAtoms("(UO2)11(CO3)6(OH)12[2-]", "C") << endl;
//	cout << ElementAtoms("(UO2)11(CO3)6(OH)12[2-]", "H") << endl;
//	cout << endl;
//	cout << ElementAtoms("Ethyne(aq)", "C") << endl;
//	cout << ElementAtoms("Ethyne(aq)", "H") << endl;
//	cout << ElementAtoms("Ethyne(aq)", "O") << endl;
//	cout << endl;
//	cout << ElementAtoms("Calcite", "Ca") << endl;
//	cout << ElementAtoms("Calcite", "C") << endl;
//	cout << ElementAtoms("Calcite", "O") << endl;
//	cout << ElementAtoms("Calcite", "H") << endl;
//}

//#include "Solubility/CO2Solubility.h"
//
//#include "Core/Multiphase.h"
//#include "Core/CanonicalReactionSystem.h"
//#include "Core/TotalEquilibriumSolver.h"
//#include "Geochemistry/GeoSystem.h"
//#include "Geochemistry/Applications/Equilibrator.h"
//#include <iostream>
//#include "Geochemistry/Assembler.h"
//
//void TotalEquilibriumSolverTest();
//
//int main()
//{
//	//DemoHaliteSolubility();
//	
//	DemoCO2Solubility();
//	
//	//TotalEquilibriumSolverTest();
//	
////	System system = AssembleSystem({"H2O", "NaCl(aq)", "Calcite"});
////	
////	Equilibrator equilibrator;
////	
////	equilibrator.AddCondition("ChargeBalanced");
////	equilibrator.AddCondition("FixedTotalMass", "H2O", 55.508);
////	equilibrator.AddCondition("FixedTotalMass", "Na[+]", 1.0);
////	equilibrator.AddCondition("FixedTotalMass", "Cl[-]", 1.0);
////	equilibrator.AddCondition("FixedTotalMass", "Ca[2+]", 2.0);
////	equilibrator.AddCondition("FixedTotalMass", "HCO3[-]", 2.0);
////	//equilibrator.AddCondition("pH", 10.060582);
////	
////	equilibrator.Equilibrate(system);
////	
////	cout << system << endl;
//	
////	System system = AssembleSystem({"H2O", "NaCl(aq)", "Calcite"}, {"Calcite"});
////	
////	system.SetComposition("Calcite", 2.0);
////	
////	Equilibrator equilibrator;
////	
////	equilibrator.AddCondition("ChargeBalanced");
////	equilibrator.AddCondition("FixedTotalMass", "H2O", 55.508);
////	equilibrator.AddCondition("FixedTotalMass", "Na[+]", 1.0);
////	equilibrator.AddCondition("FixedTotalMass", "Cl[-]", 1.0);
////	equilibrator.AddCondition("FixedTotalMass", "Ca[2+]", 1.0E-7);
////	equilibrator.AddCondition("FixedTotalMass", "HCO3[-]", 1.0E-7);
////	
////	equilibrator.Equilibrate(system);
////	
////	cout << system << endl;
////	
////	Reactor reactor(system);
////	
////	double t = 0.0, tEnd = 1000.0;
////	
////	cout << setprecision(6) << fixed;
////
////	while(t < tEnd)
////	{
////		cout << setw(15) << left << t;
////		cout << setw(15) << left << system.Get("Moles", "Calcite");
////		cout << setw(15) << left << system.Get("KineticRate", "Calcite");
////		cout << setw(15) << left << system.Get("KineticRate", "Ca[2+]");
////		cout << setw(15) << left << system.Get("pH") << endl;
////		
////		reactor.React(t, tEnd);
////	}
//	
//	
////	double tEnd = 1000.0, dt = 1.0E-3, t = 0.0;
////
////	  double T = system.GetTemperature();
////	  double P = system.GetPressure();
////	VectorXd n = system.GetComposition();
////	
////	KineticSolver kSolver(system.GetMultiphase(), system.GetEquilibriumReactions(), system.GetKineticReactions());
////	
////	kSolver.SetControlSolution(1.0E-6);
////	
////	cout << setprecision(6) << fixed;
////	
////	while(t < tEnd)
////	{
////		cout << setw(15) << left << t;
////		cout << setw(15) << left << system.Get("Moles", "Calcite");
////		cout << setw(15) << left << system.Get("KineticRate", "Calcite");
////		cout << setw(15) << left << system.Get("KineticRate", "Ca[2+]") << endl;
////		
////		kSolver.Iterate(t, tEnd, dt, T, P, n);
////		
////		system.SetComposition(n);
////	}
//	
//	//LinearInterpolatorTest();
//	
//	//LagrangeInterpolatorTest();
//	
//	//TRISpH(25.0, 1.0);
//	
//	//TRISpHTRange(Range(25.0, 150.0, 1.0), 1.0, 1.0, 0.01);
//	
//	//NoTRISpHTRange(Range(25.0, 150.0, 1.0), 1.0, 1.0E-7);
//	
//	//TestHkfModel_WaterNaCl(25.0, 1.0);
//	
//	//TestHkfModel_WaterCaCl2(25.0, 1.0);
//}
//
//void TotalEquilibriumSolverTest()
//{
//	GeoSystem system = AssembleSystem({"H2O", "NaCl(aq)"});
//	
//	Equilibrator equilibrator;
//	
//	equilibrator.AddCondition("ChargeBalanced");
//	equilibrator.AddCondition("FixedTotalMass", "H2O", 55.508);
//	equilibrator.AddCondition("FixedTotalMass", "Na[+]", 1.0);
//	equilibrator.AddCondition("FixedTotalMass", "Cl[-]", 1.0);
//	
//	equilibrator.Equilibrate(system);
//	
//	cout << system << endl;
//	
//	const Multiphase& multiphase = system.GetMultiphase();
//	
//	const CanonicalReactionSystem& reactions = system.GetEquilibriumReactions();
//	
//	BOOST_FOREACH(auto s, reactions.GetPrimarySpecies())
//		cout << s << endl;
//	
//	TotalEquilibriumSolver totalSolver(multiphase, reactions);
//	
//	VectorXd n = VectorXd::Constant(multiphase.GetSpecies().size(), 1.0E-7);
//	
//	n[system["H2O"]] = 35.508;
//	n[system["Na[+]"]] = 2.0;
//	n[system["Cl[-]"]] = 2.0;
//	
//	VectorXd uj(4);
//	
//	uj << 1.0, 55.508, -1.0E-7, 1.0;
//	
//	totalSolver.Solve(uj, 25.0, 1.0, n);
//	
//	system.SetComposition(n);
//	
//	cout << system;
//}
